library(ggplot2)
library(stargazer)
library(xtable)
library(marginaleffects)
source("functions.R")

d <- read.csv("surveydata_clean_study2.csv")

##data loading and cleaning

#delete lines with missing wave 2
d <- d[!is.na(d$disc_treatment),]

#treat immi and policy as NA for people with only one response
d$immi <- ifelse(d$immi_notna > 1, d$immi, NA)
d$disc_policy <- ifelse(d$policy_notna > 1, d$disc_policy, NA)

#imputing immigration attitudes
d$immiNA <- is.na(d$immi)
d$immi_impute <- d$immi
d$immi_impute[d$immiNA] <- mean(d$immi, na.rm=T)

#new variable for any treatment at all
d$treated <- d$disc_treatment != "control"

#new variable for treatment plus "added rhetoric" indicator
d$treat_rhet <- as.factor(paste0(d$disc_treatment, ifelse(d$rhetoric, "+rhetoric", ""), sep=""))
d$treat_rhet <- factor(d$treat_rhet, levels=levels(d$treat_rhet)[c(1,4,5,2,3)])
treatment_lvls <- levels(d$treat_rhet)

#set level ordering of treatment variable
d$disc_treatment <- factor(d$disc_treatment, levels=c("control","thematic","episodic"))


##effect of each evidenced treatment (no added rhetoric) on outcomes

#effect of each evidenced (manifest) treatment on concern (discrimination as a problem)
fit_concern <- lm(disc_problem ~ disc_treatment + immi_impute + immiNA, data=d[!d$rhetoric,])
#effect in SDs
fit_concern$coefficients / sd(d$disc_problem, na.rm=T)

#effect of each evidenced (manifest)  treatment, widespread outcome
fit_widespr <- lm(disc_widespr ~ disc_treatment + immi_impute + immiNA, data=d[!d$rhetoric,])
#effect in SDs
fit_widespr$coefficients / sd(d$disc_widespr, na.rm=T)

#Figure 2, middle panels: effect of evidenced treatments on concern and widespread
plot_estCIs(fit_concern, "evidenced_only/discr_evidenced_effects_2", "seeing discrimination as problem",
            "disc_treatment", c("thematic", "episodic"),  treatment_labels=c("audit\nstudy", "supermarket\nstory"),
            title="Study 2", ymax=.51, w=3.5)
plot_estCIs(fit_widespr, "evidenced_only/discr_evidenced_effects_widespr_2", "seeing discrimination as widespread",
            "disc_treatment", c("thematic", "episodic"),  treatment_labels=c("audit\nstudy", "supermarket\nstory"),
            title="Study 2", ymax=1, w=3.5)

#figure 3, left hand panel: effect of each evidenced treatment, policy outcome
fit_policy <- lm(disc_policy ~ disc_treatment + immi_impute + immiNA, data=d[!d$rhetoric,])
plot_estCIs(fit_policy, "evidenced_only/discr_evidenced_effects_policy_2", "favoring anti-discrimination policy",
            "disc_treatment", c("thematic", "episodic"),  treatment_labels=c("audit\nstudy", "supermarket\nstory"),
            ymax=NULL, w=3.5)

#figure 3, right hand panel: effect of each evidenced treatment, donation outcome (conditional on donating)
fit_donate <- lm(donate_Mino ~ disc_treatment + immi_impute + immiNA, data=d[!d$rhetoric & d$donate_none!=1,])
plot_estCIs(fit_donate, "evidenced_only/discr_evidenced_effects_donate_2", "donation to Mino Danmark",
            "disc_treatment", c("thematic", "episodic"),  treatment_labels=c("audit\nstudy", "supermarket\nstory"),
            ymax=NULL, w=3.5)


##additional analyses referenced in-text

#for rhetoric section: combined effect of both treatments, without rhetoric
summary(lm(disc_problem ~ treated + immi_impute + immiNA, data=d[!d$rhetoric,]) )

#for policy and donation behavior section: differential effect of the two treatments on policy
summary(lm(disc_policy ~ disc_treatment + immi_impute + immiNA, data=d[d$treated & !d$rhetoric,]))

#for policy and donation behavior section: correlates of donating
summary(lm(donate_none ~ treated, data=d))
sapply(c("disc_problem", "disc_widespr", "disc_policy", "education_int", "gender"), function(x){
        form <- as.formula(paste("donate_none ~", x))
        summary(lm(form, data=d))
      })

##analysis of comments

#merge study data with "no content" ratings, for response rate
d_treated <- d[d$treated,] #only treated respondents made comments
coding_nocontent <- read.csv("comments/surveydata_coding_nocontent.csv")
d_treated <- merge(d_treated, coding_nocontent, by="caseid")
mean(1-d_treated$no_content) #response rate, general

#merge study data and RA's coding of answers with content, for rates of emotional response
coding_affect <- read.csv("comments/surveydata_coding_affect.csv")
d_treated <- merge(d_treated, coding_affect, by="caseid", all.x = TRUE)
d_treated <- d_treated %>% mutate(emo_code_noNA = ifelse(is.na(emo_code), 0, emo_code))
#if no content, emo_code should is be 0, not NA

#rates and statistical test of emotional response (not conditional on responding)
tapply(d_treated$emo_code_noNA, d_treated$disc_treatment, mean)
summary(lm(emo_code_noNA ~ disc_treatment, d_treated))


##appendix A: representativeness

#see separate script representativeness.R


##appendix D: descriptives


##appendix E: full results behind figures

#Effect of treatments on all outcomes
stargazer(list(fit_concern, fit_widespr, fit_policy, fit_donate),
          covariate.labels=c("Audit", "Supermarket", "Immigration Att.", "Imm. Att. Missing"),
          dep.var.labels=c("Problem","Widespread","Policy","Donation"),
          omit.stat=c("f","ser"), digits=2, out="tables/fullfig345_study2.tex",
          label="tab:fullfig345_study2",
          title="Effect of the manifest-evidence treatments (without rhetoric added)
          in Study 2 on four outcomes: seeing discrimination as a problem (1-5 scale)
          and as widespread (0-10 scale), support for anti-discriminatory policies
          (1-5 scale) and donations to MinoDanmark (in DKK)")


##appendix F: full results behind figures without imputation

#Effect of treatments on all outcomes using immigration
#covariate without imputation (dropping NAs)
fit_concern_noimp <- lm(disc_problem ~ disc_treatment + immi, data=d[!d$rhetoric,])
fit_widespr_noimp <- lm(disc_widespr ~ disc_treatment + immi, data=d[!d$rhetoric,])
fit_policy_noimp <- lm(disc_policy ~ disc_treatment + immi, data=d[!d$rhetoric,])
fit_donate_noimp <- lm(donate_Mino ~ disc_treatment + immi, data=d[!d$rhetoric & d$donate_none!=1,])
fits <- list(fit_concern_noimp, fit_widespr_noimp, fit_policy_noimp, fit_donate_noimp)
stargazer(fits,
          covariate.labels=c("Audit", "Supermarket", "Immigration Att."),
          dep.var.labels=c("Problem","Widespread","Policy","Donation"),
          omit.stat=c("f","ser"), digits=2, out="tables/fullfig345_study2_noimp.tex",
          label="tab:fullfig345_study2_noimp",
          title="Effect of the manifest-evidence treatments (without rhetoric added)
          in Study 2 on four outcomes: seeing discrimination as a problem (1-5 scale)
          and as widespread (0-10 scale), support for anti-discriminatory policies
          (1-5 scale) and donations to MinoDanmark (in DKK). No imputation of
          missing immigration attitudes.")


##appendix G: weighted analyses

#Effect of each treatment on concern and widespread
fit_weight_concern <- lm(disc_problem ~ disc_treatment + immi_impute + immiNA, data=d[!d$rhetoric,], weights=weight_w2)
fit_weight_widespr <- lm(disc_widespr ~ disc_treatment + immi_impute + immiNA, data=d[!d$rhetoric,], weights=weight_w2)
fit_weight_policy <- lm(disc_policy ~ disc_treatment + immi_impute + immiNA, data=d[!d$rhetoric,], weights=weight_w2)
fit_weight_donate <- lm(donate_Mino ~ disc_treatment + immi_impute + immiNA, data=d[!d$rhetoric & d$donate_none!=1,], weights=weight_w2)
fits_weight <- list(fit_weight_concern, fit_weight_widespr, fit_weight_policy, fit_weight_donate)
stargazer(fits_weight,
          covariate.labels=c("Audit", "Supermarket", "Immigration Att.", "Imm. Att. Missing"),
          dep.var.labels=c("Problem","Widespread","Policy","Donate"), omit.stat=c("f","ser"), digits=2,
          title="Effect of all treatments in Study 2 on seeing discrimination as a problem (1--5 scale),
          perceptions of discrimination as widespread (1--10 scale), support for anti-discriminatory
          labor market policies (1--5 scale), and donations to minority charity MinoDanmark, weighting
          respondents by survey weights.",  out="tables/weighted_study2.tex")


##appendix H: heterogeneous treatment effects

#for analyses pooling Study 1 and Study see separate script discrimination_analysis_pooled.R

#analyses of additional moderators featured in Study 2 only:
outcomes <- c("disc_problem","disc_widespr")
moderators <- c("SDO", "JWB", "ethnocentrism", "RWA")
combos <- apply(expand.grid(c("Audit","Supermarket"), c("SDO", "JWB", "ethn.", "RWA")), 1, paste, collapse=", ")
ints <- lapply(outcomes, function(outcome){
  ints_permod <- lapply(moderators, function(moder){
    
    #get regression outputs
    formula <- as.formula(paste0(outcome, "~  disc_treatment*", moder))
    fit <- lm(formula, data=d[!d$rhetoric,])
    int <- tail(summary(fit)$coefficients[,c(1,2,4)], 2)
    
    #add n
    n <- nobs(fit)
    int <- cbind(int, n)
    
    #add MDE
    MDE <- int[,2]*2.8/sd(d[!d$rhetoric, outcome], na.rm=T)
    int <- cbind(int, MDE)
    
    #get CATEs at the 1st and 3rd quartile
    qs <- quantile(d[moder],  probs = c(.25, .75), na.rm=T)
    CATEs <- avg_slopes(fit, variables = "disc_treatment", by = moder)
    CATEs_at_qs <- sapply(qs, function(q){
      CATEs[CATEs[,moder]==q, "estimate"]
    })
    CATEs_at_qs <- CATEs_at_qs[c(2,1),]
    #flip order because otherwise, episode v. control is reported first
    
    #add CATEs
    int <- cbind(int, CATEs_at_qs)
    
    #shorten/adjust column names
    colnames(int) <- c("Est.", "SE", "p", "n", "MDE", "CATEq1", "CATEq3")
    
    #return results
    int
    
  })
  ints_permod <- do.call(rbind, ints_permod)
  rownames(ints_permod) <- combos
  ints_permod
})
#ints <- do.call(rbind, ints)

#add Benjamini-Hochberg adjusted p-values
pvalues <- c(ints[[1]][,3], ints[[2]][,3])
fdrs <- p.adjust(pvalues, method="BH")
ints[[1]] <- cbind(ints[[1]][,1:3], "p, BH"=fdrs[1:8], "n"=ints[[1]][,4:7])
ints[[2]] <- cbind(ints[[2]][,1:3], "p, BH"=fdrs[9:16], "n"=ints[[2]][,4:7])

#print to LateX format files
print(xtable(ints[[1]], digits=c(0,3,3,3,3,0,3,3,3)), file="tables/hetero_study2_concern.tex", only.contents = TRUE)
print(xtable(ints[[2]], digits=c(0,3,3,3,3,0,3,3,3)), file="tables/hetero_study2_widespr.tex", only.contents = TRUE)

##appendix I: effect of rhetoric on all outcomes

#effect of identity politics rhetoric on concern (seeing discrimination as a problem)
fit_rhet_concern <- lm(disc_problem ~ rhetoric + immi_impute + immiNA, data=d[d$treated,])
summary(fit_rhet_concern)

#effect of rhetoric on perceiving discrimination as widespread
fit_rhet_widespr <-lm(disc_widespr ~ rhetoric + immi_impute + immiNA, data=d[d$treated,])

#effect of rhetoric on support for anti-discrimination policy
fit_rhet_policy <-lm(disc_policy ~ rhetoric + immi_impute + immiNA, data=d[d$treated,])

#effect of rhetoric on donations to minority-related charities
fit_rhet_donate <-lm(donate_Mino ~ rhetoric + immi_impute + immiNA, data=d[d$treated & d$donate_none!=1,])

#in a table
fitlist <- list(fit_rhet_concern, fit_rhet_widespr, fit_rhet_policy, fit_rhet_donate)
outcomes <- c("disc_problem","disc_widespr","disc_policy","donate_Mino")
rhet_coefs <- t(sapply(1:4, function(i){
  
  #get coefs and n
  fit <- fitlist[[i]]
  coefs <- summary(fit)$coefficients["rhetoricTRUE",c(1,2,4)]
  n <- nobs(fit)
  
  #add MDE
  outcome <- outcomes[i]
  if(outcome=="donate_Mino"){
    MDE <- coefs[2]*2.8/sd(d[d$treated & d$donate_none!=1, outcome], na.rm=T)
  } else {
    MDE <- coefs[2]*2.8/sd(d[d$treated, outcome], na.rm=T)
  }
  
  #return
  c(coefs, n=n, MDE=MDE)
  }))
rownames(rhet_coefs) <- c("Problem","Widespread","Anti-discr. policy","Donation")
colnames(rhet_coefs) <- c("Est.", "SE", "p", "n", "MDE")

#print to Latex
print(
  xtable(rhet_coefs, digits=c(0,3,3,3,0,3),
         caption="Effect of adding identitarian rhetoric to the evidenced discrimination treatments on four outcomes (Study 2).",
         label = "tab:rhetoric"),
  file="tables/rhetoric.tex"
)
